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_  _--  ABSTRACT 

The  computer  program  for  the  analysis  of  linearly  elastic  plane- 
stress  or  plane-strain  problems  devised  by  Felippa  in  his  work  on 
"Refined  Finite  Element  Analysis  of  Linear  and  Nonlinear  Two-dimensional 
Structures"  has  been  modified  to  include  the  use  of  initial  displacement 
boundary  conditions.   In  addition  the  original  IBM  7094  computer  dependent 
program  has  been  adapted  for  use  on  the  IBM  360/65  computer.   In  both 
progrsms  the  FORTRAN  IV  language  has  been  used. 

Problems  involving  "Poor  fit"  displacement  boundary  conditions  and 
refined  mesh  analysis  using  coarse  mesh  analysis  input  displacements, 
which  could  not  have  been  done  with  the  original  program,  are  now  pos- 
sible with  the  modified  version  presented  herein. 
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Chapter  I     Introduction 


In  this  chapter  reasons  are  given  for  the  use  of  computer  structural 
analysis.   The  purpose  of  the  presentation  is  explained,  and  some  expected 
program  uses  are  listed.   The  objectives  and  extent  of  work  completed  by 
the  author  are  presented. 


1. 1  Reasons  for  computer  analysis 

The  engineer  of  the  present  day  is  faced  with  solving  structural 
problems  of  great  complexity.   Classical  mathematics,  despite  its  ever 
increasing  sophistication,  is  only  capable  of  solving  severely  idealized 
situations  while  at  the  same  time  placing  a  burden  on  skilled  manpower 
which  could  be  better  used  for  design  and  development  processes. 

Fortunately,  the  simultaneous  development  of  the  digital  computer  and 
new  general  methods  of  solution  for  problems  in  Continuum  Mechanics  has 
come  to  the  rescue  in  many  areas  of  investigation.   The  modern  digital  com- 
puter coupled  with  the  finite  element  method  is  revolutionizing  the  ap- 
proach to  the  process  of  analysis.   Expensive  and  time  consuming  experi- 
mental models,  now  often  used  in  the  design  of  important  structures,  are 
rapidly  becoming  displaced  by  more  economical  computation. 

1.2  Purpose  of  presentation 

This  presentation  offers  a  computer  program  which  takes  advantage  of 
the  modern  computer's  high  computational  speed  and  the  versatility  of  the 
finite  element  method.   The  program  is  for  the  general  analysis  of  arbif 
trary  plane-stress  or  plane-strain  structural  mechanics  problems  in  linear 
elasticity.   Sufficient  background  material  on  both  the  finite  element 
method  and  the  program  structure  is  given  so  that  those  not  conversant 
with  these  matters  may  use  the  program  effectively.   In  addition,  informa- 
tion on  the  practical  aspects  of  the  program  alogrithm  for  the  direct  stif- 
fness procedures,  boundary  condition  application,  and  solution  of  large 
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systems  of  linear  equations  are  given  for  those  interested  in  making 
modifications  or  additions  to  the  program. 

1.3  Expected  uses 

Two-dimentional  elastic  analysis  by  the  finite  element  method  is 
well  suited  for  the  development  of  large  scale,  production  usage  computer 
programs.   Expected  areas  of  use  for  the  program  presented  include: 

(1)  Practical  working  design  situations,  ranging  from  one-time 
preliminary  analyses  to  extensive  in-depth  studies. 

(2)  Design  data  generation  for  stress  concentration  factors  or 
design  tables  for  specific  structural  shapes. 

(3)  Augmentation  of  other  analysis  methods.  Localized  area 
analysis.  Utilization  of  the  thermal  stress  and  body  force 
features.  Gross  analysis  to  determine  instrument  or  sensor 
placement . 

1.4  Scope  and  objectives 

The  original  version  of  the  program  was  written  by  Dr.  Carlos  A. 
Felippa  [3]*  during  the  course  of  his  Ph.D.  studies  at  the  University  of 
California.   The  original  program  was  dependent  on  the  IBM  7094  computer 
and  associated  hardware  then  at  the  Berkeley  Computer  Center. 

The  extent  of  the  author's  objectives  and  work  on  the  program  has 
been:  (1)  reprogramming  modifications  necessary  for  the  adaptation  of  the 
program  to  the  faster  IBM  360/65  computer  at  the  Naval  Postgraduate  School 
Computer  Center,  (2)  substituting  the  use  of  the  more  expedient  random  ac- 
cess disk  unit  for  the  original  tape  storage  methods  which  were  used  for 


* 
Numbers  in  brackets  refer  to  the  list  of  references  on  page  5^ 


12 


portions  of  the  external  storage  requirements,  (3)  restructuring  and 
modifying  the  program  to  allow  usage  of  initial  displacement  boundary 
conditions,  thus  allowing  the  program  to  accept  all  physical  boundary 
conditions  for  planar  problems,  (4)  creation  of  a  "User's  Manual"  to 
facilitate  use  of  the  program  by  persons  not  familiar  with  the  finite 
element  method  or  production  programming  techniques.   This  presentation 
is  orientated  toward  this  final  objective. 
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Chapter  II     General  Program  Features 


In  this  brief  chapter  the  general  program  capabilities  and  features 
are  given  as  background  information. 


2.1  Capabilities 

The  program  in  its  current  form  represents  a  large  scale,  high  speed 
general  computational  processor  for  the  stress  analysis  of  plane  elastic 
bodies.   It  is  one  of  the  few  such  programs  available  in  the  open  litera- 
ture. 

The  program  is  capable  of  analyzing  any  arbitrary  plane  structural 
shape,  singly  or  multiply  connected.   The  body  thickness  may  vary,  to  a 
degree.   All  physical  boundary  conditions  are  acceptable.   Structural 
loading  types  may  include  surface  forces  (concentrated  and  distributed), 
body  forces,  and  thermal  loading  effects.   Any  linearly  elastic,  homogen- 
ous, isotropic  material  may  be  analysed.   Up  to  six  different  materials 
are  allowable  for  composite  material  structures. 

2.2  Analytic  results 

The  method  of  analysis  is  the  finite  element  procedure  [2],  [16], 
with  the  displacement  solution  method  and  direct  stiffness  matrix  genera- 
tion [12],  [13].   Displacement  compatible  finite  elements  are  used.   This 
type  of  element  has  the  virtue  in  theory  that  if  finer  and  finer  mesh  sub- 
divisions are  used  and  there  is  no  numerical  error  or  round  off,  conver- 
gence to  the  exact  solution  is  assured  [11], 

In  the  analysis  the  program  computes  in-plane  deflections  and  stresses 
at  selected  sites  on  the  body,  resulting  from  in-plane  loading.   Stress 
components  Cf  $       ff  ,    and   7"   are  computed,  as  well  as  the  principal 
stresses  and  directions  at  every  nodal  point.   Stress  contour  graphs  of 
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the  structure,  or  a  magnified  subportion,  are  generated,  wherein  the 
structure  outline  and  constant  stress  level  contour  lines  are  printed, 


15 


Chapter  III     Fundamentals  of  the  Finite  Element  Method 


In  this  chapter,  the  fundamentals  of  the  finite  element  method  are 
given  in  order  to  acquaint  the  user  with  the  procedures  and  the  nomen- 
clature needed  to  understand  the  program  structure  and  fully  utilize  its 
potentials. 


3. 1  Basic  idealization 

The  basic  concept  of  the  plane  finite  element  method  is  that  any 
continuous  two-dimensional  body  may  be  separated  by  imaginary  lines  into 
a  finite  number  of  individual  elements,  i.e.,  an  assemblage  of  smaller 
individual  plates.   Thus  an  actual  continuous  structure  is  replaced  by 
an  assemblage  of  discrete  structural  elements  (Fig.  3.1). 
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Figure  3.1   The  Finite  Element  Idealization 
The  elements  are  assumed  to  be  interconnected  only  at  a  discrete 
number  of  joints  or  nodal  points  through  which  forces  are  transmitted, 
but  the  structural  elements  are  such  that  the  connectivity  of  the  struc- 
tural system  is  preserved  along  all  the  common  boundries  of  adjacent 
elements. 

With  these  stipulations:  (1)  continuum  approximation  by  an  assem- 
blage of  elements  and  (2)  element  displacement  field  restrictions,  we 
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have,  in  effect,  converted  a  continuous  body  into  a  body  composed  of  many 
structural  elements  where  each  element  can  deform  in  only  a  certain  number 
of  predictable  shapes. 

A  given  point  in  the  continuum  has  two  independent  degrees  of  freedom 
of  displacement.   The  corresponding  point  in  the  idealized  structure  still 
has  two  degrees  of  freedom,  but  the  displacement  field  is  now  restricted 
to  the  one  selected  and  imposed  on  the  element.   This  reduction  is  basic 
to  any  solution.   In  the  finite  element  method  the  reduction  is  made  in  a 
physical  manner  using  structural  elements  that  can  retain  many  of  the 
physical  properties  and  behavior  characteristics  of  the  original  material. 
In  other  methods  of  analysis  the  reduction  is  often  made  in  a  mathematical 
procedure,  in  which  many  of  the  desirable  properties  and  characteristics 
are  lost  through  inability  to  make  adequate  or  tractable  mathematical 
models. 
3.2  Displacement  fields 

The  constraints  we  have  placed  on  the  element  take  the  form  of  dis- 
placement functions,  or  fields,  for  which  the  overall  element  displace- 
ment shapes,  and  resulting  strains  and  stresses,  are  functionally  re- 
lated to  the  nodal  point  displacements  (generalized  coordinates).   The 
success  of  the  element  in  duplicating  the  performance  of  the  actual  con- 
tinuous body  is  critically  dependent  on  the  choice  of  the  displacement 
function.   Much  effort  has  been  given  to  the  determination  of  successful 
displacement  functions  and  the  fundamental  rules  governing  their  genera- 
tion. For  flat  elements,  general  procedures  have  been  developed  [3],  how- 
ever as  yet  there  is  no  such  general  method  for  the  construction  of  dis- 
placement fields  for  a  curved  element  [1], 
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3.3  Convergence  requirements 

To  guarantee  that  a  finite  element  solution  will  converge  to  the 
true  solution,  under  conditions  of  increasingly  finer  element  division 
(finer  mesh),  the  planar  displacement  function  must  be  able  to  satisfy 
three  basic  conditions: 

(1)  Continuity  condition:  interelement  displacement  compatability 
must  be  satisfied.   For  plane  cases  this  implies  that  adjacent 
element  edges  must  remain  together  under  element  deformation. 

(2)  Completeness  conditions:  (a)  rigid  body  modes  -  overall 
rigid  body  displacement  of  the  element  must  not  result  in 
element  straining;  (b)  constant  strain  states  -  nodal  displace- 
ments that  are  indicative  of  constant  strain  must  result  in  con- 
stant strain  conditions  in  an  element. 

(3)  Invariance  -  element  stiffness  properties  derived  through 
the  use  of  displacement  functions  must  remain  the  same  for  all 
coordinate  systems  which  are  used  for  their  derivation. 

3.4  Element  equilibrium  equations  and  stiffnesses 

Knowing  element  displacement  modes  as  functions  of  nodal  displacements 
we  can  now  write  the  equilibrium  equations  for  the  element.   The  relations 
between  nodal  displacements  and  nodal  forces  (generalized  forces)  are  ex- 
pressed in  compact  form  by  the  element  stiffness  matrix.   The  resulting 
element  equilibrium  equations  take  the  form: 


1 


(3.1) 


Where, 


|  F>   =  Nodal  point  force  vector  for  an  element,  (2N  x  1) 

^  e 
1(1   =  Nodal  point  displacement  vector  for  an  element,  (2N  x  1) 
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[k]   ■  Element  stiffness  matrix,  symmetric,  (2N  x  2N) . 
(N  is  the  number  of  nodes  per  element) 

The  element  stiffness  matrix  is  a  function  of  the  geometric  and  physical 
properties  of  the  element. 

The  derivation  of  element  stiffnesses  can  be  approached  using  the 
principle  of  virtual  displacement  [13],  [16],   Here,  virtual  nodal  dis- 
placements in  a  local  frame  of  reference  are  imposed  on  the  element.   The 
external  and  internal  work  done  by  the  various  forces  and  stresses  during 
displacement  are  equated.   The  resulting  equation  can  be  reduced  to  the 
general  form: 


where, 


F    = 


[*fWM*H 


(3.2) 


[B]  =  Strain-Displacement  relationships  from 
elasticity  considerations. 

[D]  =  Stress-Strain  relationships  from  elementary 
theory. 

The  element  stiffness  is  developed  by  performing  the  indicated  matrix 
integration. 

An  alternate  method  of  deriving  the  element  stiffness  is  an  energy 
approach  [3],  [51 •   Here  the  total  potential  energy  of  the  element-load 
system  is  minimized,  an  outgrowth  of  the  application  of  the  governing 
variational  principle.   This   is,  in  fact,  a  form  of  the  Ritz  technique 
applied  to  the  network  of  finite  elements.   This  method  is  not  as  direct, 
nor  as  physically  interpretable  as  the  previous  approach.  However  it  is 
a  more  powerful  technique  and  leads  naturally  to  the  formulation  of  the 
consistent  mass  matrix,  consistent  stiffness  matrix  and  consistent  load 
vectors.   Elements  having  sophisticated  geometries  and  higher  order  dis- 
placement fields  can  be  handled  practically  only  with  the  energy  method. 


19 


Having  the  element  stiffness  matrix  and  the  two  auxiliary  matrices 
[B]  and  [D]  we  can  develop  the  full  stress  analysis  for  a  single  element 
when  the  nodal  displacements  resulting  from  the  loading  are  known. 

3.5  Complete  stiffness  matrix 

To  permit  combining  the  individual  stiffnesses  each  must  first  be 
transformed  from  its  local  coordinate  system  to  a  common  (global)  system 
for  the  entire  node-element  mesh.   This  is  done  with  a  conventional  tensor 
transformation  relating  the  two  coordinate  systems.   The  formulation  of 
the  complete  stiffness  matrix  for  the  discretized  structure  can  now  be 
accomplished.   In  the  direct  stiffness  method  this  is  achieved  by  the 
direct  addition  of  the  element  stiffnesses  at  all  the  element  interface 
nodal  points.   This  is  a  technique  idealy  suited  for  computer  operation 
(see  appendix  2). 

3.6  System  equilibrium  equations 

The  complete  stiffness  matrix  [K]  is  the  matrix  of  equilibrium  equa- 
tions for  the  total  idealized  structure.   The  set  of  overall  equilibrium 
equations  is  of  the  same  form  as  equation  (3.1),  viz.: 


K 


(3.3) 


The  nodal  point  force  vector  (load  vector)  contains  the  loads  to  which  the 
body  is  subjected.   The  loads  are  in  the  form  of  concentrated  nodal  point 
forces  directed  along  the  global  coordinates  axes.   Where  the  body  loading 
is  in  the  form  of  distributed  forces,  body  force  loads  or  thermal  loads, 
equivalent  concentrated  forces  are  applied  at  the  appropriate  nodal  points. 

The  nodal  point  displacement  vector  (displacement  vector)  is  the  un- 
known in  equation  (3.3).   This  is  a  direct  result  of  the  formulation 
developed  and  is  responsible  for  the  name,  "Displacement  Solution  Method." 
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Another  general  method  of  finite  analysis  is  available,  where  a  stress 
field  is  used  [5],  this  procedure  is  called  the  "Equilibrium  Method." 

The  overall  equilibrium  equations  must  be  modified  by  the  boundary 
conditions  of  the  problem.   This  involves  elimination  of  stiffness  contri- 
butions in  [K]  for  constrained  nodes.   In  the  case  of  initial  displace- 
ment boundary  conditions  the  load  vector  and  displacement  vector  must  be 
adjusted  to  reflect  equivalent  nodal  forces  and  the  initial  displacement 
value,  respectively,  that  result  from  displacing  a  node  a  known  amount. 

The  final  set  of  equations  (3.3)  can  attain  rather  large  proportions. 
The  order  of  the  load  and  displacement  vectors  and  the  size  of  the  sym- 
metric overall  stiffness  matrix  are  2  x  P,  where  P  is  the  total  number  of 
nodal  points  in  the  mesh.   Fortunately  the  stiffness  matrix  is  sparsely 
populated  and  banded,  which  permits  specialized  solution  techniques  on  a 
computer.   Actually,  the  displacement  solution  method  and  the  direct  stif- 
fness procedures  are  used  so  that  these  stiffness  matrix  properties  may 
be  exploited  (see  appendix  2). 

The  final  equations  (3.3)  are  solved  for  the  nodal  point  displace- 
ments.  With  the  known  displacements,  the  matrices  [B]  and  [D]  may  be 
applied  for  each  element  and  the  overall  nodal  stresses  evaluated.   The 
stress  levels  at  nodes  common  to  more  than  one  element  are  usually  averaged 
to  obtain  a  representative  value.   The  end  result  is  the  displacements 
and  stress  levels  at  each  node  of  the  idealized  structure. 
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Chapter  IV     Program  Information  and  Structure 


This  chapter  offers  information  on  the  program  and  gives  its 
structure,  functional  routines  and  general  method  of  operation. 


4.1  Program  identification 

PSELST  -  Plane  Stress  Elastic  Analysis  using  Linear  Strain  Triangles. 
Programmed:   Carlos  A.  Felippa,  June  1966 
Revised:     John  P.  Malone,  July  1968 

4.2  Purpose 

The  purpose  of  the  program  is  to  provide  a  high  speed,  production 
use  computational  solution  of  general  plane-stress  or  plane-strain  static, 
linear  elastic  problems  using  linear  strain  triangles  in  the  finite  ele- 
ment method.   Surface  loads,  body  forces  and  thermal  effects  may  be  con- 
sidered. 

4.3  Programming  information 

The  program  is  written  in  FORTRAN  IV  language  [6]  for  the  IBM  OS/360 
Model  67  computer.   An  overlay  structure  [7],  [8],  under  control  of  the 
system  linkage  editor,  is  utilized  to  conserve  main  core  storage.   The 
overlay  segments  are  arranged  so  that  only  active  subroutines  and  atten- 
dant internal  storage  requirements  are  in  main  core  at  a  time, 

4.4  External  storage 

Fortran  logical  units  1,2,3,7,8,  and  9  are  used  for  temporary  stor- 
age.  Logical  unit  7  is  defined  as  a  data  set  on  one  IBM  2311  random  ac- 
cess disk  storage  unit.   The  remaining  logical  units  are  defined  as 
separate  data  set  blocks  on  a  second  2311  disk  unit.   The  input  card 
reader,  output  printer,  and  output  card  punch  are  defined  as  units  5,6, 
and  4  respectively. 
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4.5  Basic  finite  element  mesh  units 

The  basic  mesh  element  is  a  quadrilateral  composed  of  four  6-nodal 
point  linear  strain  triangles  (LST) ,  the  center  point  being  the  centroid, 
(Fig.  4.1).   Internal  points  9  to  13  are  eliminated  by  matrix  condensation 
[3],  [9]  thereby  reducing  the  number  of  degrees  of  freedom  from  26  to  16. 

Single  6-nodal  point  triangles  may  also  be  specified  to  facilitate 
fitting  of  certain  shapes. 


KEY: 


Quadrilateral 
(4  corner  nodes,  8  external  nodes) 
(#9  ■  centroid) 


Corner  Nodes 

Midside  Nodes 

Removed  during 
condensation 


Single  Triangle 
(3  corner  nodes,  6  external  nodes) 


Figure  4.1  Basic  Finite  Element  Mesh  Shapes 
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4.6  Capacity 

The  mesh  input  is  subject  to  the  following  limitations  for  an  IBM 
computer  with  256  K  bytes  of  main  core  storage; 

Max.  number  of  elements  (Quad,  and/or  triangle)  -  -  350 
Max.  number  of  external  (total)  nodal  points  -  -  -  1050 
Max.  number  of  restrained  components  -------   250 

Max.  difference  of  nodal  point  numbers  for 

the  same  element  -----------------   79 

These  limits  are  dictated  by  storage  requirements  in  the  computation 
of  the  overall  stiffness  matrix,  and  not  by  the  equation  solver. 

It  is  estimated  that  with  512  K  bytes  of  main  core  storage  the  pro- 
gram could  be  modified  to  accept  1700  external  nodes  and  a  half-band- 
width of  250.   This  would  result  in  approximately  600  usable  elements. 
Solution  times  for  full  scale  problems  with  these  modifications  may  be 
prohibitive  on  current  computers. 

The  most  commonly  encountered  limitation  of  the  current  program  ver- 
sion is  the  maximum  nodal  point  number  difference.   Next  in  order  are 
maximum  external  nodes  (total  nodes),  maximum  number  of  elements,  and 
maximum  number  of  restrained  components. 

The  maximum  number  of  degrees  of  freedom  is  1050  x  2 (external  nodes) 
+  350  x  10(condensed  and  recovered  nodes)  ■  5600.   The  condensation  pro- 
cedure reduces  the  number  of  equations  to  2100,  yet  retains  the  versatil- 
ity of  5600  degrees  of  freedom. 

The  maximum  half-bandwidth  (2  x  Max.  node  number  difference  +  2) ,  a 
measure  of  the  width  of  the  banded  area  of  the  overall  stiffness  matrix, 
is  160. 
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4.7   Program  structure 

The  overlay  root-segment  structure  is  shown  in  Figure  4.2,  where 
each  rectangle  represents  a  subroutine.   The  subroutine  functions  are: 


MAIN 

RDDISK- 
WRDISK 


SETUP 

STQUAD 
STLST6 
LDINPT 

THERLD 

FORMK 

SOLVE 

BIGSOL 

STRESS 

TRISTR 

CNTPLT 


remains  in  core  and  controls  calling  sequence  during  execution; 

remains  in  core  and  provides  high-speed  deposit  and  retrieval 

of  data  on  a  random  access  disk  storage  unit; 

inputs,  prints  and  checks  mesh  data,  and  evaluates  element 

stiffnesses; 

assembles  and  condenses  quadrilateral  stiffness; 

computes  stiffness  of  a  six  nodal  point  triangle; 

inputs  load  cases  and  reduces  surface,  body,  and  thermal  loads 

to  equivalent  forces  on  element  external  nodal  points; 

computes  initial  thermal  forces  for  a  single  triangle; 

assembles  the  complete  stiffness  matrix; 

obtains  nodal  point  displacements  from  BIGSOL; 

solves  large  capacity  banded  matrix  problems; 

evaluates  and  prints  element  and  nodal  stresses; 

computes  stresses  for  a  single  triangle; 

produces  printer  plots  of  stress  contour  lines. 


The  subroutine  flow  chart  is  presented  in  Fig.  4.3. 
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-  START  - 


V 


SETUP 


All  transfers 
indicated  by 

are  controlled  by 
MAIN,  and  require 
loading  and 
unloading  of 
overlay  segments. 


-+y< 


yes. 


LDINPT 


^st  I 
Load  | 
Case  | 


FORMK 


SOLVE 


STRESS 


I 
I 


CNTPLT 


I 


ADDED 
LOAD 


no| 


N 


PROBLEM 


no 


^Sj) 


STQUAD 


I 


STLST6 


THERLD 


Added 

Load 

Cases 

I 


BIGSOL 


\ 


TRISTR 


S 
\ 

I 

No 

Plots 
I 
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Figure  4.3   Subroutine  Flow  Chart 
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4.8   Programming  techniques 

Three  programming  techniques  that  are  contained  in  the  original 
version  of  the  program  are  noteworthy: 

(1)  RDDISK/WRDISK  Subroutine  -  used  as  a  data  transfer  device  to 
external  storage.   Benefits  are:  simplified  data  address  accounting 
in  subroutines,  fast  data  transfer,  access  to  large  storage  areas. 

(2)  Matrix  Compression  (condensation)  -  upper  half-band  of  the 
overall  stiffness  matrix  is  stored,  column-wise,  by  blocks  for 
ease  of  manipulation  and  to  reduce  storage  requirements. 

(3)  Use  of  one  dimensional  array  -  to  minimize  address  calculation 
time. 

Techniques  (2)  and  (3)  are  described  more  completely  in  appendix  2. 
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Chapter  V     User  Information 


In  this  chapter  specific  features  of  the  program  are  explored, 
Techniques  are  given  to  assist  in  the  effective  use  of  the  program. 


5.1  Typical  applications 

Representative  applications  of  specific  problem  types  are  illustrated 
in  Fig.  5.1.   Case  (a)  is  a  wing  frame  where  gust  loading  is  simulated  by 
a  static  load  situation.   Thermal  loading  effects  at  the  frame  tip  might 
also  be  considered.   Case  (b)  is  a  cross-section  of  a  length  of  pressure 
piping.  The  circular  perimeters  have  been  distorted  during  a  bending 
process.   The  stress  levels  and  patterns  at  the  thin  wall  section  are 
desired.   Case  (c)  is  a  comparison  study  of  two  welds.   The  left  fillet 
welds  are  cheaper  to  produce  in  some  given  process,  but  the  right  side 
full-penetration  welds  are  expected  to  be  significantly  stronger.  A 
cost  effectiveness  study  requires  information  comparing  the  two  welds. 
Case  (d)  is  a  stress  concentration  factor  study.  Data  for  design  tables 
is  desired  to  guide  design  engineers  confronted  with  this  configuration. 
Loading  situations  of  tension,  compression  and  in-plane  bending  might  be 
considered.   Case  (e)  is  a  section  of  a  prefabricated  tunnel  proposed  for 
underwater  installation.   The  tunnel  is  partially  evacuated  as  part  of  a 
high  speed  transportation  scheme.   The  optimum  tunnel  cross  section  is 
desired,  using  a  minimum  material  criterion.   Case  (f)  is  a  hypothetical 
plane-strain  problem  illustrating  the  versatility  of  the  finite  element 
technique  used  in  conjunction  with  a  high  speed  computer. 

5.2  General  program  input  and  output 

The  versatility  of  the  program  is  obtained  by  accepting  a  penalty  of 
requiring  a  large  amount  of  input  data.   The  data  is  in  the  form  of  stand- 
ard punched  cards.   The  number  of  input  cards  will  be  between  25  and  500, 
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dependent  on  the  number  of  elements  used  and  the  extent  of  loading. 

Input  information  is  divided  in  two  groups,  structure  data  and  load- 
ing data.   The  general  content  of  input  data  follows: 

1.  Structure  data. 

a.  Accounting  totals  -  number  of  elements,  nodes,  etc. 

b.  Mesh  configuration. 

1.  Node-to-element  identification. 

2.  Corner  nodes,  X-Y  coordinates. 

c.  Material  physical  constants. 

d.  Boundary  conditions 

2.  Loading  data. 

a.  Accounting  totals  -  number  and  types  of  loads. 

b.  Load  magnitudes  and  locations. 

Once  a  problem  has  been  set  up  with  the  structure  data  input,  repeated 
load  cases  may  be  run.   Consecutive  problems  may  also  be  run. 

Program  output  is  divided  into  three  parts:  (1)  echo  check  of  input 
data,  (2)  displacements  and  stresses,  (3)  contour  maps.   The  general  con- 
tent of  each  part  follows: 

1.  Echo  check. 

a.  Input  structure  data  printed  for  each  problem. 

b.  Input  loading  data  printed  for  each  load  case. 

2.  Displacements  and  stresses. 

a.  Final  load  vector. 

b.  Nodal  point  displacements  printed  (punch  option). 

c.  Computed  element  nodal  stresses  (optional). 

d.  Averaged  nodal  stresses,  (J"  ,  rf  ,  T    »  CT      >  C7".  , 
&  '  ^x'  ^y    xy    max    min 

7"   ,  and  maximum  principle  stress  orientation  (punch 
max 

option) . 
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3.   Contour  graphs 

a.  A  contour  graph  for  each  averaged  stress  computed 
(optional) . 
Details  of  input,  output  and  program  options  will  be  covered  in  subsequent 
sections. 
5.3   Precautions 

This  section  is  presented  here  to  prevent  users  from  rushing  to  a 
card  punch  machine  with  this  manual  in  their  hand. 

The  author  has  found  that  the  following  "rules"  must  not  be  violated 
when  using  the  program. 

(1)  Know  your  problem. 

(2)  Be  patient. 

(3)  Be  prepared. 

An  example  of  the  first  rule  involves  how  much  of  the  structure  to 
analyze.   Many  structures  have  1-fold  symmetry,  less  frequently  2-fold. 
If  the  loading  is  symmetric  about  the  same  axes,  partial  body  analysis 
may  be  utilized.   This  avoids  duplication  of  data  and  allows  more  ef- 
fective use  of  the  elements  available. 

The  rule  on  patience  should  be  applied  during  mesh  construction, 
nodal  point  numbering,  and  card  punching.   The  program  has  many  error 
printout  exits  that  determine  and  illuminate  input  mistakes,  but  these 
exits  also  terminate  program  execution. 

Being  prepared  means  having  a  plan  to  make  effective  use  of  the 
tremendous  amount  of  information  that  the  program  produces.   A  cursory 
glance  at  the  sigma-max  graph  to  determine  if  a  material  is  suitable  may 
fail  to  make  adequate  use  of  the  wealth  of  output  information  available. 
Other  persons  may  need  information  from  your  analysis  or  some  unusual 
condition  may  exist  that  warrents  further  study. 
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5.4  Mesh  construction 

An  effective  mesh  is  one  that  represents  a  good  approximation  to  the 
true  structure.  How  effective  a  mesh  is  may  not  be  known  until  one  run 
has  been  completed  and  stress  patterns  are  reviewed. 

An  efficient  mesh  is  one  that  is  formed  on  the  basis  of  the  above  re- 
quirements, but  also  allows  computation  to  proceed  as  rapidly  as  possible. 
This  section  deals  with  the  construction  of  efficient  meshes. 

The  mesh  idealization  of  the  body  may  contain  1  to  350  basic  mesh 
units  (Fig.  4.1).  Due  to  its  superior  stiffness  properties,  use  of  the 
quadrilateral  is  preferred.   The  use  of  triangles  should  be  limited  to 
curve  fitting  difficulties  and  when  changing  the  fineness  (or  coarseness) 
of  the  mesh;  this  is  termed  grading. 

The  nodes  of  the  basic  mesh  unit  are  called  corner  points  or  midside 
points.   Mesh  units  may  only  be  joined  with  correspondance  between  corner 
and  midside  points.  Two  adjoining  units  share  one  midside  point  and  two 
corner  points.   In  a  mesh,  corner  points  will  join  three  (usually  four)  or 
more  elements  when  internal  to  the  body  (Fig.  5.2c,  points  a  and  b).   Mid- 
side  points  are  always  common  to  only  two  elements  when  within  the  body 
interior.   Each  element  and  node  must  be  numbered  for  identification  pur- 
poses.  Element  size,  shape  and  location  information  is  determined  by 
listing  the  X-Y  (global)  coordinates  of  all  corner  points.  Midside  point 
coordinates  and  element  local  coordinate  systems  are  internally  generated 
in  the  program.  Nodal  numbering  techniques  and  coordinate  system  usage 
are  covered  in  the  following  sections. 

Three  factors  should  be  considered  in  constructing  the  mesh: 

(1)  Expected  or  known  stress  gradients. 

(2)  Mesh  symmetry. 

(3)  Nodal  point  numbering. 
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Stress  gradients  should  dictate  the  degree  of  fineness  of  the  mesh. 
In  areas  of  high  stress  gradients  the  mesh  should  be  refined.   Portions 
where  stress  levels  are  nearly  constant  can  be  represented  with  large 
elements  in  a  coarse  mesh.   In  many  problems  users  intuitatively  produce 
meshes  that  account  for  stress  gradients.   However,  when  complex  shapes, 
boundary  conditions  and  loads  are  involved,  a  preliminary  coarse  mesh 
analysis  should  be  conducted.   The  indication  of  high  stress  gradients  is 
not  always  available  from  the  contour  graphs,  since  they  are  based  on 
averaged  values  of  common  nodal  stresses.   A  true  indication  of  the  need 
for  element  refinement  is  when  element  stress  levels  at  shared  nodes 
differ  considerably. 

Many  techniques  for  grading  the  mesh  during  initial  construction,  or 
when  refining  is  required,  are  available.   A  comparison  of  two  grading 
methods  is  shown  in  Figure  5.4. 

Mesh  symmetry  with  respect  to  some  line  or  curve  in  the  body  should 
be  maintained  where  possible.   In  complex  shapes  overall  symmetry  usually 
cannot  be  achieved.   In  these  cases  local  symmetric  areas  should  be  con- 
structed.  Figure  5.2  illustrates  symmetric  mesh  construction. 

In  a  completed  mesh  the  one  element  that  has  the  largest  difference 
between  its  nodal  point  numbers  influences  heavily  the  computation  time 
required  for  the  problem.   The  maximum  difference  in  node  numbers  is  an 
indication  of  the  amount  of  storage  required  to  house  and  solve  the  set 
of  equilibrium  equations.   In  solving  the  equations  the  data  in  storage  is 
transferred,  manipulated,  and  restored  many  times.  The  larger  the  data 
set,  the  longer  the  corresponding  computation  time.   The  mesh  can  be  con- 
structed so  that  it  lends  itself  to  optimum  nodal  point  numbering  tech- 
nique, which  in  turn  will  reduce  the  largest  nodal  point  number  difference 
and  the  computation  time. 
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(a) 
Unsymmetric  (poor)  Mesh 


Improved  Symmetric  Mesh 
(reduced  computation  time) 


KEY: 

^Xf      Corner  Nodal  Point 


Midside  Nodal  Point 


(b) 
Symmetric  Mesh 


Figure  5.2   Symmetry  in  Mesh  Construction 
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In  general,  the  least  nodal  number  difference  can  be  achieved  by  plac- 
ing all  mesh  corner  nodes  on  lines  or  smooth  curves  within  the  body.   Figure 
5.2c  illustrates  this  procedure. 

It  is  suggested  that  the  physical  process  of  mesh  construction  be  ac- 
compolished  in  the  following  manner.  (1)  Sketch  initial  mesh  ideas  on  a 
scratch  pad  until  the  final  version  is  formed.  (2)  Check  the  limitations 
on  maximum  numbers  of  elements,  nodes,  and  constraints  and  the  maximum 
nodal  number  difference.  (3)  Transfer  the  mesh  layout  to  a  large  sheet 
where  sufficient  space  is  available  to  identify  corner  and  midside  nodes 
by  number  and  to  exhibit  the  X-Y  coordinates  of  corner  nodal  points.   Large 
(2'x3')  Ozalid  prints  from  a  master  drawing  of  faint  %  inch  square  grid 
with  a  superimposed  quadrant  of  concentric  arcs  have  been  found  convenient 
for  this  last  step. 

Meshes  can  be  produced  where  all  the  elements  are  of  equal  size  and 
have  the  same  orientation.   In  this  case  all  elements  have  the  same  stif- 
fness matrix  and  the  program  need  compute  only  one.   This  is  a  time  saving 
feature  that  can  be  used  for  regular  bodies  or  preliminary  gross  analyses. 
Frequently  an  auxiliary  computer  program  can  be  written  to  generate  all 
mesh  input  data  cards  in  these  cases. 
5.5  Nodal  point  and  element  limbering 

The  program  requires  that  each  element  and  each  nodal  point  be  numbered 
for  identification  purposes.   The  numbering  of  elements  and  nodes  is  in- 
dependent, but  increasing  both  counts  in  the  same  pattern  is  convenient. 
Consecutive  numbering  is  not  a  requirement,  but  again,  it  is  convenient. 
If  consecutive  numbering  is  not  used  (occurs  when  combining  two  previously 
independent  meshes)  an  ordered  listing  of  the  input  sequence  should  be 
maintained,  since  output  is  identified  by  consecutive  numbering. 
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The  method  of  nodal  point  numbering  is  important  for  two  reasons. 

(1)  To  keep  within  the  limit  on  maximum  nodal  point  difference 
for  a  single  element. 

(2)  Efficient  numbering  can  significantly  reduce  computation 
time. 

The  general  rule  for  efficient  numbering  is  to  number  in  a  direction 
where  the  least  number  of  nodes  lie  in  "line"  (not  necessarily  the  small- 
est figure  dimension)  and  repeat  for  the  next  "parallel  line."  The 
technique  is  illustrated  in  figures  5.3  and  5.4c. 

22 
17 
14 

1  o 


2   i, 


Figure  5.3   Nodal  Numbering  Technique 
Efficient  numbering  within  graded  areas  is  often  difficult,  Figure  5.4 

compares  two  grading  techniques  and  three  numbering  methods.   Figure  5.4(c) 

is  preferred  on  both  scores. 

As  discussed  previously,  the  one  element  in  the  mesh  with  the  greatest 
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node  number  difference  controls  the  solution  time.   If  it  is  unavoidable 
that  a  large  nodal  number  difference  exists  in  some  area  of  the  mesh, 
nothing  is  to  be  gained  by  trying  to  reduce  lower  difference  values  else- 
where in  the  mesh. 

5.6  Coordinates 

The  global  coordinate  system  for  the  program  is  the  planar  cartesian 
coordinate  system.   The  X  and  Y  coordinates  of  each  corner  point  are  re- 
quired input  data.   This  coordinate  array  data  set  is  usually  the  largest 
single  block  of  input  data.  Two  hints  are  offered  that  can  reduce  the 
time  to  prepare  these  data  cards.  (1)  Reduce  the  significant  figures  of 
the  coordinate  by  slight  nodal  point  movement.   This  procedure  may  disrupt 
mesh  symmetry,  but  the  effect  will  be  small.   (2)  The  coordinate  origin 
can  be  placed  internal  to  the  figure.   Mesh  symmetry  may  be  utilized  to 
produce  identical  coordinate  values,  except  for  sign. 

The  coordinate  system  is  compatible  with  structure  overlapping.   This 
situation  has  not  been  explored  other  than  to  demonstrate  its  dreadful  in- 
fluence on  the  contour  graphs,  which  should  not  be  printed  in  these  cases. 
Overlapping  can  also  be  produced  as  a  result  of  deformation  and  the  pro- 
gram does  not  sense  the  condition. 

5.7  Units 

A  consistent  set  of  units  must  be  used  between  the  linear  measure  and 
surface  traction  load  values.  Where  the  inch  is  used  as  the  unit  of  bound- 
ary and  thickness  measure,  distributed  forces  and  shears  must  be  per  square 
inch  (eg.  PSI) .   The  modulus  of  elasticity  and  specific  weight  input  con- 
stants must  reflect  the  same  units.   The  most  useful  set  of  units  for 
mechanical  design  has  been  the  inch  and  KSI  (i.e.  Kip  per  square  inch);  for 
large  structural  problems  the  foot  and  KSF. 
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5.8  Materials 


The  program  is  written  for  homogenous,  isotropic,  linear  elastic 
materials.   The  required  input  data  are  the  Elastic  Modulus,  E;  Poisson's 
Ratio,  1/     ;  Specific  Weight,  0   J  and  Coefficient  of  Thermal  Expansion,  <X  . 
For  plane-strain  problems  reduced  (primed)  values  must  be  used: 

The  ability  to  utilize  multiple  materials  is  one  of  the  advantages 
of  the  finite  element  procedure.   In  the  program  up  to  six  different 
materials  may  be  specified.   True  difference  between  materials  may  be 
used,  or  different  property  values  of  the  same  material.   Temperature 
influence  on  the  elastic  properties  may  be  approximated  by  using  "differ- 
ent" materials  when  in  the  thermal  gradient  of  a  single  material  object. 
Body  force  loading,  other  than  the  program  standard  of  1-G,  can  be  simu- 
lated by  applying  suitable  multipliers  to  the  value  of  the  specific  weight 
prior  to  input.   This  technique  is  useful  in  simulating  dynamic  inertia 
loading  by  a  static  equivalence. 

When  multi-material  problems  are  analyzed,  the  element  stress  print 
option  should  be  specified.   The  standard  print  of  averaged  nodal  point 
stresses  is  not  adequate  at  material  interfaces.   The  stress  contour  graphs 
should  be  used  only  with  judgement  in  these  cases,  since  they  do  not  re- 
flect true  stress  discontinuities  at  interfaces. 
5.9  Boundary  conditions 

Boundary  conditions  are  applied  at  the  appropriate  nodal  points 
(corner  or  midside)  of  the  structure.  A  boundary  constraint  consists  of 
fixing  a  node  in  one  direction.  A  "fixed"  boundary  point  will  require  two 
boundary  constraints.   The  program  will  accept  up  to  250  constraints. 
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Constraints  are  not  restricted  to  nodes  on  the  physical  boundary;  how- 
ever, this  is  the  usual  situation. 

Initial  values  of  boundary  displacement  may  also  be  used  as  boundary 
conditions.   This  feature  has  many  applications.   One  is  the  simulation 
of  "poor  fit"  conditions  in  a  part  or  structure.  A  second  is  illustrated 
in  Figure  5.5,  where  actual  gear  tooth  root  area  deflections  were  deter- 
mined; then  used  as  input  initial  displacements,  thus  freeing  elements  for 
use  in  a  more  refined  mesh  of  the  structure  of  interest. 

In  program  input  the  boundary  condition  data  is  designated  for  the 
appropriate  nodal  point  with  a  boundary  condition  "TAG."  A  list  of  TAG 
values  and  appropriate  boundary  conditions  follow: 
TAG  VALUE  BOUNDARY  CONDITION 

TAG  =0  If  the  point  is  fixed  in  both  directions. 

or 

Initial  X  and  Y  displacements  are  specified  and 
the  point  is  fixed  in  both  directions. 
TAG  =1  If  the  point  is  fixed  in  the  X-direction  and 

free  to  move  in  the  Y-direction. 
or 

Initial  X-displacement  is  specified  but  the  point 
is  free  to  move  in  the  Y-direction. 
TAG  =2  If  the  point  is  free  to  move  along  a  line  forming 

an  angle  (D     with  the  positive  X-axis, 
or 

Initial  Y-displacement  is  specified  but  the  point 
is  free  to  move  in  the  X-direction. 
The  boundary  conditions  are  illustrated  in  Figure  5.6. 
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Figure  5.5  Use  of  Initial  Displacement  Boundary  Conditions 
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Figure  5.6   Boundary  Conditions 
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A  single  initial  displacement  can  only  be  specified  in  the  X-direc- 
tion  (TAG=1)  or  the  Y-direction  (TAG=2,  (f)=   0.0°).   Initial  displacements 
for  points  on  a  slope  roller  are  not  available  in  this  version  of  the  pro- 
gram.  There  is  an  equivalance  between  TAG=1  and  TAG=2,<j9  =  90.0°.   The 
first  method  requires  less  internal  computation. 

In  multi-load  case  problems  any  initial  displacements  specified  will 
apply  only  for  the  first  load  case.   Subsequent  load  cases  will  have  the 
same  constraint  conditions  for  initially  displaced  nodes,  but  with  a  0.0 
value  of  initial  displacement. 

Boundary  constraints  are  a  fundamental  characteristic  of  a  structure 
and  are  known  immediately  upon  the  definition  of  a  problem.  When  the  ad- 
vantages of  2-fold  or  4-fold  symmetry  are  utilized,  with  partial  body 
analysis,  the  constraints  needed  are  not  always  obvious.   Care  should  be 
exercised  in  these  cases  to  ensure  duplication  of  the  actual  problem 
situation.   When  a  complete  structure  is  analyzed  that  is  shape  and  load 
symmetric,  auxiliary  boundary  constraints  are  often  helpful  to  improve  the 
solution.   Frequently  symmetric  centerline  nodes  can  be  placed  on  rollers. 
This  will  avoid  slight  skewing  displacements  which  result  from  round-off 
error  in  the  stiffness  matrix  generation. 
5.10  Loading 

Three  loading  types  are  available  with  the  program:  (1)  surface  loads, 
(2)  body  forces,  and  (3)  thermal  loading.   The  general  character  of  the 
input  loading  data  is  listed  below; 

Loading  type  Input  characteristics 

Concentrated  force  Magnitude  of  X  and  Y  components  specified. 

Loaded  node  identified. 

Distributed  normal  force*     Magnitude  and  sense  of  normal  pressure. 

Element  side  identified. 

Distributed  shear  force*      Magnitude  and  sense  of  surface  shear. 

Element  side  identified. 
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Gravity  (body  force)  load     Indicated  true  (T)  if  applicable. 

Acts  in  (-Y)  direction. 

Thermal  load  Indicate  temperature  increment  at 

element  corner  nodes. 

*  -  side  variation  assumed  parabolic  based  on  values  specified 
at  the  nodes  of  the  element  side. 

The  sign  convention  for  concentrated  forces  is  that  of  the  global 
coordinate  system.  Distributed  forces  follow  a  sign  convention  based  on 
the  element  shape.   Normal  outward  traction  and  counterclockwise  shears 
are  positive  (Fig.  6.2).   Gravity  loading  acts  in  the  (-Y)  direction. 

Distributed  forces  and  shears  acting  on  an  element  side  are  assumed 
to  have  a  parabolic  variation  based  on  the  input  values  at  the  corner 
nodes  and  midside  node  involved.  The  resulting  total  element  side  load 
is  represented  by  equivalent  nodal  point  concentrated  forces  that  are 
generated  internally  in  the  program.  When  the  gravity  load  option  is  used, 
the  program  generates  nodal  forces  equivalent  to  body  forces  to  apply  at 
all  the  mesh  nodes.   In  thermal  loading  the  program  computes  the  forces 
that  would  be  required  to  reestablish  the  size  of  the  thermally  expanded 
elements  and  applies  them  at  the  appropriate  nodes.  These  computations 
are  internal  to  the  program,  with  no  external  indication  other  than  the 
resulting  displacements  and  stresses. 
5.11  Contour  graphs 

The  stress  contour  plots  generated  by  the  program  are  very  useful. 
The  plots  are  produced  on  a  conventional  printer,  thus  eliminating  a  re- 
quirement for  specialized  plotting  equipment  not  available  at  many  computer 
installations.   The  loss  in  definition  due  to  printer  plotting  is  compen- 
sated for  by  the  short  time  required  to  produce  the  plots  (about  9  sec). 

Six  plots  are  available,  each  shows  the  body  outline  and  contour  lines 
of  constant  stress.   The  lines  divide  the  full  range  of  stresses  encountered 
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into  ten  equal  increments.   Plots  of  the  following  stresses  are  available; 

CTX,  (Ty»        7^y»   ^tax»  ^in'   ^max* 

The  plotting  subroutine  allows  up  to  50  elements  to  be  eliminated 

from  the  area  considered  in  determining  the  stress  range  for  the  contour 

lines.   This  feature  allows  (1)  the  plotting  of  only  a  portion  of  the 

structure,  (2)  elimination  of  high  stress  gradient  areas  that  bunch  the 

countour  lines,  (3)  a  combination  of  the  two. 

5.12  Error  exits 

Several  error  conditions  caused  by  either  wrong  input  data  or  exceeded 
program  capacity  are  checked  in  subroutine  SETUP.   Error  messages  printed 
are  self-explanatory  and  may  be  complemented  by  examination  of  the  input 
data  printout.   The  program  does  not  stop  until  all  error  conditions  have 
been  tested.   If  another  problem  follows  after  an  error  detection,  the 
program  searches  for  the  next  START  indicator  card  at  which  time  execu- 
.  tion  continues  on  the  new  problem. 

5. 13  Displacement  solution  iterations 

Double  precision  iterations  for  improvement  of  the  displacement  solu- 
tion are  optionally  available  in  the  program.   One  or  two  iterations 
should  always  be  used  for  large  ill-conditioned  problems.   Each  iteration 
is  essentially  re-solving  the  displacement  equations  to  eliminate  round- 
off error,  and  will  take  approximately  90%  of  the  original  problem  solu- 
tion time. 

5. 14  Card  punching  of  displacements  and  stresses 

A  program  option  is  the  punching  of  nodal  displacements  and  averaged 
nodal  point  stresses.   These  card  data  sets  may  be  used  as  input  to  plot- 
ting equipment  for  the  production  of  more  sophisticated  plots  than  pro- 
duced by  the  program. 
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5.15  Timing 

Sufficient  problems  have  not  been  run  with  the  program  to  develop  an 
accurate  empirical  relationship  giving  overall  execution  time.   An  esti- 
mator currently  in  use  is: 

60  nodes,  half -bandwidth  <  60 

1  minute  execution  time  for  each:       50  nodes,  60  -^  half -bandwidth  <  120 
(no  iterations) 

40  nodes,  120  ^  half -bandwidth 

For  small  problems  (number  of  elements  <  40)  execution  times  pre- 
dicted will  be  low  since  the  subroutine  SETUP  requires  a  minimum  of  50 
sec.  regardless  of  problem  size.   Timing  estimates  apply  to  an  IBM  0S/360 
computer. 

5.16  Accuracy 

To  illustrate  the  accuracy  of  this  version,  a  2:1  symmetric  plate 
under  uniform  in-plane  load  on  top  (Fig.  5.7)  and  constrained  by  dia- 
phragms at  the  ends  (X-disp.  only)  was  analyzed  by  subdividing  one-half 
of  the  plate  into  meshes  of  (n  x  n)  quadrilaterals  -  square  in  this  case. 
There  is  no  known  analytic  solution,  but  this  problem  has  been  used  ex- 
tensively in  comparison  tests  of  different  types  of  finite  elements  be- 
cause of  the  diversity  of  stress  conditions  and  the  ease  of  preparation. 
Precision  is  not  sought  in  the  analysis,  since  a  "true"  solution  is 
guaranteed  if  the  mesh  is  sufficiently  refined.   Rapidity  of  convergence 
is  the  desirable  characteristic. 

A  comparison  of  the  typical  values  reproduced  in  Figure  5.7  shows 
that  the  8x8  mesh  provides  almost  5  decimal  digits  for  the  displace- 
ments, 4  for  C7"  and  3  for  7"   •  Actually,  the  program  has  capacity  for 
an  18  x  18  mesh  (5500  degrees  of  freedom)  if  necessary.   The  consistency  of 
the  stress  values  is  reflected  by  the  fact  that  the  maximum  discrepancy  over 
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contributing  elements  nodal  points  did  not  exceed  0.004  for  C     and 

rf    .  and  0.006  for   7"    in  the  case  of  the  8x8  mesh. 
u  y9  xy 
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Chapter  VI     Input  Data  Preparation 


This  chapter  is  written  as  a  self  contained  unit  giving  all  specific 
information  for  the  preparation  of  input  data  cards.   When  facility  with 
the  use  of  the  program  has  been  attained,  this  chapter  may  be  duplicated 
and  used  as  a  condensed  program  user's  manual. 


FLOATING  POINT  (F)  FORMATS  ARE  COMPATIBLE  WITH  ASSIGNED  EXPONENT  (E)  FORMATS 
RIGHT  JUSTIFY  INTEGER  AND  EXPONENT  NUMBERS  IN  THEIR  ASSIGNED  FIELDS 
6. 1  Structure  data 

(a)  Start  Card  (A8) :  With  the  word  START  punched  in  cols.  1-5. 
This  card  must  precede  the  input  data  deck  of  any  problem. 

(b)  Title  Card  (20A4) :   Alphameric  information  in  cols.  1-80  to 
identify  the  output. 

(c)  Control  Card  (814,  6L2) : 

Variable 
Columns  Name  Meaning 

1-4  NUMEL  Number  of  elements  (^  350); 

5-8  NUMCP  Number  of  corner  points; 

9-12  NUMNP  Number  of  total  nodal  points  (>$  1050); 

13-16  NUMBC  Number  of  restrained  nodal  points; 

17-20  NUMPB  Number  of  plot  defining  boundary  points, 

see  item  (e) ; 

21-24  NLOAD  Number  of   load   cases;   will  be    set=l   if 

left  blank; 

25-28  NMAT  Number  of  different  materials  (-$6); 

will  be  set=l  if  left  blank; 

29-32  MAXIT  Maximum  number  of  residual  iterations 

in  the  displacement  solution;  punch  a 
1  or  2  for  large,  ill-conditioned  prob- 
lems. 

The  next  five  fields  are  for  logical  flags;  if  a  T  is  punched,  in  any 

assigned  column,  the  indicated  option  takes  place.   A  blank  of  F  implies 

FALSE. 
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33-34  Tl  All  quadrilaterals  have  the  same  stiffness 

matrix  (see  note  1) ; 

35-36  T2  Punching  of  mesh  nodal  point  coordinates 

and  displacements  (14,  2F8.3,  2E14.5); 

37-38  T3  Punching  of  averaged  (X  %     C  ,   and  i 

at  mesh  nodal  points  and  quadrilateral 
centroids  (14,  3E18.6); 

39-40  T4  Print  of  element  stresses  (see  note  2); 

41-42  T5  Another  complete  problem  follows. 

Notes 

(1)  All  quadrilaterals  have  the  same  stiffness  if  they  can  be  super- 
imposed by  a  translation. 

(2)  Element  stresses  should  be  printed  in  problems  involving  several 
material  types,  since  averaged  stresses  and  their  plots  do  not 
display  actual  interface  discontinuities. 

(d)  Material  Property  Table  (14,  E10.3,  2F10.3,  E10.4):  One  card  per 
material  type  (total  NMAT  cards) : 

Cols.     1-4      Material  type  number; 
5-14     Elastic  modulus; 
15-24     Poisson's  ratio; 
25-34     Specific  weight; 
35-44     Coefficient  of  thermal  expansion. 
For  plane-strain,  reduced  values  must  be  used; 

£'=7i%        v'=r^w      <*'=*(n-y) 

(e)  Plotted  figure  boundary  outline  array  (2014) :  For  the  plotting 
of  stress  graphs,  NUMPB  number  of  corner  nodal  points  which,  when  con- 
nected outline  the  structure,  or  structure  subregion,  as  a  series  of 
straight  lines  must  be  punched,  in  cyclic  order,  20  node  numbers  per  card. 
The  starting  corner  node  and  direction  of  travel  around  the  structure,  or 
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subregion,  boundary  are  arbitrary.   Holes  in  multiply  connected  bodies 
cannot  be  outlined  separately.   In  subregion  plots,  the  elements  (^  50) 
outside  the  plotted  figure  boundary  (skipped  elements)  must  be  listed  in 
item  (1)  under  Loading. 

(f)  Element  nodal  point  array  (1014,  F10.3):   One  card  per  element. 
(total  NUMEL  cards). 

Cols.     1-4      Element  number; 

5-36     Nodal  point  numbers: 

(I)  for  quadrilateral:   element  nodal  points  in 
counterclockwise  order  I-J-K-L-M-0-P  (Fig.  6.1a). 

t 

The  starting  corner  is  arbitrary,  except  when 
equal  stiffnesses  are  implied  (i.e.,  T1»T  in  Control 
Card).   In  this  case  the  starting  corners  must  be 
in  the  same  location  for  all  elements. 

(II)  for  single  triangles:   punch  nodes  I-J-K-L-M-N, 
(Fig.  6.1b),  leave  cols.  29-36  blank. 

37-40     Element  material  type,  will  be  set=l  if  left  blank. 

41-50     Element  thickness,  will  be  set=] .0  if  left  blank. 
Note :    If  a  quadrilateral  is  not  convex  (not  recommended)  the  entrant 
corner  must  be  either  J  or  K. 

(g)  Corner  point  coordinate  array  (14,  2F8.3):   One  card  per  corner 
point  (total  NUMCP  cards). 

Cols.     1-4      Corner  nodal  point  number; 
5-12     X-coordinate 
13-20     Y-coordinate 
(h)  Boundary  condition  array  (214,  2E15.3):   One  card  per  restrained 
nodal  point  (total  NUMBC  cards). 
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Cols.     1-4      Nodal  point  number 

5-8      Tag  =  0  if  point  fixed  in  both  directions  or  initial 
displacement  is  specified  in  both  the  X  and  Y 
direction  then  point  is  fixed. 
=  1  if  point  is  fixed  in  the  X-direction  and  free 
in  the  Y-direction  or  initial  displacement  is 
specified  in  the  X-direction  and  point  is  free 
in  the  Y-direction. 
=  2  if  point  is  free  to  move  along  a  line  forming 
angle  (D     with  the  X-axis  and  fixed  in  a  direc- 
tion normal  to  that  line;  or  if  initial  dis- 
placement is  specified  in  the  Y-direction  and 
the  point  is  free  in  the  X-direction. 
9-23     Angle  in  degrees,  positive  counterclockwise  from  X-axis 
for  type  TAG=2  boundary  condition. 

Initial  X-displacement  value  for  TAG=0  boundary  condi- 
tion (for  both  cases  will  be  set=0.0  if  left  blank) 
24-38     Initial  displacement  boundary  condition  value.  Y- 

displacement  for  TAG=0  or  2.  X-displacement  for  TAG=1. 
(will  be  set  =  0.0  if  left  blank) 
6.2  Loading  data 

Each  load  case  must  be  specified  by  a  data  deck  initiated  by  a  LOADING 
card;  this  package  follows  the  structure  data  deck.  A  load  deck  consists 
of  the  following  cards; 

(i)  Loading  Card  (A8) :   With  the  word  LOADING  punched  in  Cols.  1-7. 
(j)  Title  Card  (20A4) :   Alphameric  information  in  cols.  1-80  to 
identify  the  load  case. 
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Figure  6.1   Element  Nodal  Point  Identification 
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Figure  6.2   Convention  for  Element  Side  Loading 
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(k)  Control  card   (314, L2) 
Variable 
Columns     Name  Meaning 

1-4      NPLD  Number  of  nodal  points  with  concentrated 

forces; 
5-8      NELD  Number  of  element  sides  loaded  with  dis- 

tributed forces; 
9-12     NTLD  Number  of  elements  undergoing  thermal 

increments; 
13-14     DENS  Logical  flag  for  gravity  loading:   If  a  T 

(implies  .TRUE.)  is  punched,  gravity  forces 
acting  along  the  (-Y)  direction  are  con- 
sidered. 
(1)  Stress  Contour  Graph  Indicator  Card  (714):  A  positive  integer 
punched  in  any  of  the  first  six  fields  will  cause  a  stress  contour  plot 
to  be  produced,  on  the  printer,  for  the  indicated  stress  component. 
Columns  Graph 

1-4  Sigma  x 

5-8  Sigma  y 

9-12  Tau  xy 

13-16  Sigma  max 

17-20  Sigma  min 

21-24  Max  shear 

The  last  field  (cols.  25-28)  indicates  the  number  NSK  -^  50  of  ele- 
ments to  be  skipped  from  the  plots.   If  NSK  >  0,  additional  cards  must 
follow,  specifying  the  skipped  element  numbers  (2014).   Element  skipping 
may  be  used  for  two  different  purposes: 

(1)  to  eliminate  small  regions  of  high  stress  gradients  which 
cannot  be  accurately  described  by  a  printer  plot; 
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(2)  to  plot  a  portion  of  the  structure,  which  is  then  amplified. 

In  this  case  the  plotted  figure  outline  array  (item  e)  must 

specify  the  outline  of  the  subregion. 

(m)  Nodal  Point  Forces  (14,  2F8.3):  One  card  per  nodal  point  loaded 

with  a  concentrated  force  (no  cards  if  NPLD  ■  0). 

Cols.      1-4      Nodal  point  number; 

5-12     X-load; 

13-20     Y-load. 

(n)  Element  Side  Loads  (214,  6F8.3):  One  card  per  element  side  under 

surface  traction  (no  cards  if  NELD  =  0).   The  convention  for  positive 

traction  and  shear  is  indicated  in  Figure  6.2.   The  side  variation  is 

assumed  to  be  parabolic  and  specified  by  the  values  at  points  a,  b  and  c 

(in  counterclockwise  sense).   For  instance,  for  side  2  of  a  quadrilateral: 

a  =  corner  point  J,  b  =  midside  point  N,  c  =  corner  point  K. 

Cols.     1-4      Element  number; 

5-8      Side  number  (see  Fig.  6.2); 

9-16     Normal  traction  at  a  (T  )  ; 

a 

17-24  Normal   traction  at  b    (T.); 

b 

25-32     Normal  traction  at  c  (T  ) : 

c 

33-40     Surface  shear  at  a; 
41-48     Surface  shear  at  b; 
49-56     Surface  shear  at  c. 
These  values  must  be  specified  per  unit  of  length  of  the  figure  and 
per  unit  of  thickness,  (eg„,  PSI) . 

(o)  Thermal  Increments  (14,  4F10.3):   One  card  per  element  under- 
going temperature  changes  (no  cards  if  NTLD  =  0): 


56 


Cols.     1-4      Element  number; 

5-14     Temperature  variation  at  corner  I; 
15-24     Temperature  variation  at  corner  J; 
25-34     Temperature  variation  at  corner  K; 
35-44     Temperature  variation  at  corner  L  (leave  blank  for 
a  triangle). 
Note;   the  thermal  increment  at  the  centroid  of  quadrilaterals  is  assumed 
to  be  the  mean  of  the  corner  values,  and  a  linear  variation  assumed  over 
each  subtriangle. 
6.3  New  problem 

The  input  of  a  new  problem  must  follow  the  last  load  deck  for  the 
previous  one.   For  safety,  any  number  of  blank  cards  may  be  inserted 
before  the  START  card. 
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Recommendations 

The  most  hearty  recommendation  is  to  use  the  program.   It  represents 
a  powerful  engineering  tool. 

Recommendations  for  augmentation  of  the  method  are: 

(1)  Inclusion  of  an  optional  bilinear  element  stiffness  subroutine 
to  allow  analysis  of  the  ever  increasing  group  of  materials  that 
have  such  characteristics. 

(2)  Modification  to  accept  linearly  varying  element  thickness. 

(3)  A  general  orthotropic  element  stiffness  subroutine. 
Recommendations  for  augmentation  of  the  program  are: 

(1)  A  mesh  generation  package  configured  for  the  NPS  Computer 
Center  IBM  Optical  Display  Unit. 

(2)  A  contour  graph  plotting  package  for  use  on  commercial  X-Y 
plotters,  and  compatible  with  the  program  punched  card  output  of 
displacements  and  stresses. 

(3)  Conversion  of  tape-disk  external  storage  statements  in  the 
program  to  direct  access  statements. 
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APPENDIX  1     Sample  Problem  -  "Flounder"  plate  in  tension 


This  appendix  presents  a  sample  problem  analyzed  with  the  program. 
The  structure  and  loading  situation  are  somewhat  hypothetical  to  allow 
presentation  of  a  variety  of  input  data  and  computed  results.  Actual 
input  data  decks  and  computer  output  are  illustrated. 


A  "Flounder"  plate  in  the  Naval  service  is  any  roughly  triangular 
plate,  that  in  its  lifetime  may  inadvertently  find  itself  on  the  ocean 
floor.  The  example  problem  presented  in  this  appendix  is  the  analysis 
of  one  type  of  flounder  plate. 

The  function  of  the  problem  plate  is  that  of  an  attachment  member. 
The  plate  is  welded  to  a  "Trolley-Block"  apparatus  which  rides,  on  cables, 
between  ships  at  sea  conducting  underway  transfer  of  supplies.   The 
flounder  plate  acts  as  the  attachment  point  for  a  small  fixture,  which  in 
turn  holds  nets  or  boxes  containing  the  transferred  material.   The  problem 
is  the  preliminary  coarse  mesh  analysis  of  the  plate. 

Sketches  of  the  complete  trolley-block  assembly  and  the  loaded 
flounder  plate  are  presented  in  Figure  Al.l.   The  figure  also  illustrates 
the  idealized  form  of  the  plate;  nodal  point  numbers;  element  numbers; 
nodal  coordinates;  boundary  conditions,  and  input  loading. 

The  program  input  data  deck  is  illustrated  in  Figure  A1.12.   Data 
values  that  are  overlined  in  the  figure  represent  input  for  which  optional 
default  values  are  available.   The  size  of  the  data  fields  for  input  data 
are  also  illustrated. 

The  computed  output  data  is  presented  in  Figure  A1.3,  where  signi- 
ficant items  are  annotated  on  the  figure. 
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APPENDIX  2     Practical  Aspects  of  the  Computer  Alogrlthm 

This  appendix  presents  practical  aspects  of  the  computer  alogrithm 
for  the  direct  stiffness  procedures,  equation  solver,  and  the  applica- 
tion of  boundary  conditions. 


The  direct  stiffness  procedure  is  the  process  whereby  the  complete 
structure  stiffness  matrix  is  formed  from  the  individual  element  stiffness 
matrices.   The  process  is  conceptually  simple,  since  it  requires  only  the 
systematic  addition  of  the  stiffnesses  of  all  the  elements  in  the  system. 
Difficulties  arise  when  the  concept  is  transposed  to  a  computer  process 
where  storage  space  is  at  a  premium.   A  firm  understanding  of  the  methods 
of  assembly  and  storage  of  the  complete  stiffness  matrix  within  the  pro- 
gram is  of  importance  for  those  desiring  to  modify  the  existing  program. 
These  methods  control  the  programming  techniques  for  application  of  bound- 
ary conditions,  and  the  solution  of  the  equilibrium  equations. 

The  superposition  of  each  element  stiffness,  in  the  formation  of  the 
complete  stiffness  matrix,  is  accomplished  by  adding  its  individual  terms 
into  the  complete  stiffness  matrix  according  to  the  nodal  point  number 
of  the  elements.   This  is  illustrated  in  Figure  A2.1. 

The  program  could  theoretically  construct  the  complete  stiffness 
matrix  in  the  manner  illustrated  in  Fig.  A2.1;  a  two-dimensional  array 
could  be  defined  and  simple  indexing  methods  utilized  to  superpose  the 
element  stiffnesses.   The  storage  requirements  for  this  straightforward 
approach  would  be  a  square  array  with  side  dimension  of  (2  x  Maximum 
Number  of  Nodes).   The  resulting  array  (2100  x  2100)  would  occupy  the 
total  main  core  storage  of  70  computers  of  the  size  used  for  this  pro- 
gram.  The  techniques  used  to  alleviate  the  storage  problem  utilizes 
three  characteristics  of  the  complete  stiffness  matrix. 

The  complete  stiffness  matrix  is 
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Sfiffness   Coefficient 


Element     Stiffness     Matrix 
(  16  x   16  ) 


FIG.    A2.I        DIRECT      STIFFNESS      PROCEDURE 
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(1)  Symmetric 

(2)  Banded 

(3)  Sparsely  populated 

These  properties  are  illustrated  in  Figure  A2.2. 


Sparseness 
(Zero  values) 


Figure  A2.2   Complete  Stiffness  Matrix  Characteristics 
Symmetry  permits  a  reduction  of  approximately  one-half  in  the  storage 
required,  since  only  a  triangular  half  of  the  matrix  is  required  to 
retain  significant  data.   The  banded  and  sparseness  properties  allow 
further  savings,  since  only  data  between  the  diagonal  and  the  band  limits 
is  non-zero  and  required  in  computation.   The  cross-hatched  portion  of 
Fig.  A2.2  shows  the  actual  quantity  of  storage  required.   The  storage  space 
problem  is  only  partially  solved,  since  space  approximately  equal  to  an 
array  size  dimensioned  (%  Bandwidth  x  Diagonal  Length)  is  still  required,, 
In  the  program  the  maximum  half -bandwidth  is  160  and  the  maximum  diagonal 
length  (equals  a  side  dimension)  is  2100.   The  resulting  (160  x  2100)  ar- 
ray is  small  enough  to  be  stored  on  an  external  disk  storage  unit,  but 
cannot  be  housed  completely  in  main  core  (5  would  be  required)  where  high 
speed  calculations,  independent  of  internal/external  storage  transfer 
time,  can  be  performed. 
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Considerations  of  the  equation  solving  method  and  transfer  time 
reduction  dictate  the  final  techniques  used  in  the  storage  and  manipula- 
tion of  the  data  contained  in  the  complete  stiffness  matrix.   These  two 
considerations  will  be  discussed  before  proceeding. 

The  equation  solver  used  in  the  program  employs  Cholesky's  algorithm 
[4],  [10] •   The  fundamental  concept  of  the  solution  process  is  to  perform 
operations  on  the  complete  stiffness  matrix  to  reduce  it  to  triangular 
form,  after  which  the  unknown  displacements  are  found  by  back  substitu- 
tion procedures.   The  limited  coupling  characteristic  of  the  set  of 
equations  allows  the  reduction  in  form  to  be  accomplished  in  a  manner 
where  a  subgrouping  of  equations  can  be  reduced  independently  in  a  pro- 
cess termed  triangular  decomposition.  The  decomposition  process  does 
not  modify  the  load  vector  in  the  equilibrium  equations,  and  once  com- 
pleted, the  decomposed  form  of  the  complete  stiffness  matrix  may  be 
stored  for  possible  use  with  additional  load  cases.  This  accounts  for 
the  greatly  reduced  running  time  for  additional  load  cases.   The  sub- 
groups of  equations  are  handled  as  individual  blocks  of  equations. 
Figure  A2.3  shows  a  block  structure  used  for  decomposing  the  complete 
stiffness  matrix.  The  figure  also  illustrates  two  points  concerning 
block  structure;  (1)  element  stiffness  contributions  may  enter  two  dif- 
ferent blocks,  and  (2)  in  order  to  maintain  equal  block  size,  the  first 
must  contain  a  zero  value  area. 

The  transfer  time  between  external  and  internal  core  storage  is  a 
function  of  the  data  address  determination  time  and  the  physical  time 
required  to  transfer  the  data.  Where  many  individual  units  or  small 
groups  of  data  are  being  moved  from  within  a  larger  data  set  (the  situa- 
tion in  our  case),  the  address  computation  time  is  usually  the  largest 
contributor  to  the  overall  time.  To  reduce  address  time  to  a  minimum, 
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zero  value  area 


Block 
Numbers 
Circled 


A  single  element  may 
contribute  to  two 
— ,    blocks 


Figure  A2.3  Block  Structure  used  for  decomposing  [K] 
the  use  of  singly  dimensioned  arrays,  with  their  attendent  short  address 
calculations,  is  required. 

Returning  to  the  complete  stiffness  matrix,  with  its  maximum  data 
storage  requirement  of  a  (160  x  2100)  array,  the  considerations  of  the 
preceeding  paragraphs  are  applied.   The  final  storage  technique  for  the 
resulting  matrix  data  calls  for  data  in  block  form  with  each  block  stored 
individually  on  an  external  disk  unit.   Each  block  is  transferred,  as 
required,  to  main  core  storage  for  the  formation,  modification,  or  solu- 
tion procedures.   Each  block  is  stored  in  a  one-dimensional  array  form. 
The  overall  storage  method  is  illustrated  in  Figure  A2.4,  where  an  orig- 
inal complete  stiffness  matrix  is  shown,  along  with  a  virtual  array 
depicting  the  actual  data  that  requires  storage,  and  the  final  one- 
dimensional  array  form  of  one  of  the  stored  blocks.   Block  divisions  are 
indicated  in  the  figure.   The  rows  of  the  virtual  array,  and  the  order  of 
data  in  the  one -dimensional  array,  are  the  columns  (read  bottom  to  top) 
of  the  upper  half-bandwidth  of  the  complete  stiffness  matrix.   The  first 
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column  of  the  virtual  array  is  the  main  diagonal  of  the  complete  matrix. 
The  original  rows  become  diagonals  in  the  virtual  array. 

The  application  of  boundary  conditions  to  the  structure  equilibrium 
equations  involves  adjustment  of  the  complete  stiffness  matrix  to  re- 
flect the  constrained  nodes  or  initially  displaced  and  constrained  nodes. 
In  the  case  of  constrained  nodes,  the  condition  is  accounted  for  by  strik- 
ing out  the  rows  and  columns  corresponding  to  the  degrees  of  freedom 
associated  with  the  boundary  condition  and  replacing  the  corresponding 
diagonal  terms  with  a  non-zero  value.   If  initial  displacements  are  speci- 
fied the  corresponding  column  of  the  matrix  must  be  multiplied  by  the 
initial  displacement  value,  and  the  resulting  vector  subtracted  from  the 
load  vector  before  the  above  constraint  procedures  are  carried  out.   In 
addition,  the  known  displacement  value  must  be  inserted  into  the  displace- 
ment vector.   The  boundary  condition  application  procedure  is  illustrated 
in  Figure  A2.4,  where  a  5-element  bar,  fixed  at  the  ends  (nodes  3  and  28); 
with  an  initial  positive  Y-displacement  at  node  13,  is  loaded  with  a 
negative  Y-direction  concentrated  force,  P,  at  node  16. 

The  carry  over  of  boundary  conditions  into  the  virtual  storage  array 
and  one  block  of  the  final  storage  form,  is  also  illustrated  in  Figure 
A2.4.   Since  the  stiffness  data  is  actually  stored  in  a  singly  dimensioned 
array,  special  indexing  and  address  accounting  procedures  (see  subroutine) 
FORMK)  must  be  used  to  allow  complete  row  and  column  boundary  condition 
procedures  to  be  applied  in  situations  where  two  data  blocks  are  involved; 
each  block  being  available  for  manipulation  only  when  in  main  core. 
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Appendix  3 


Program    (PSELST)   Listing 
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